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ABSTRACT 


Particle-based hydrodynamics models offer distinct advantages over Eulerian and 
Lagrangian hydrocodes in particular shock physics applications. Particle models are designed 
to avoid the mesh distortion and state variable diffusion problems which can hinder the 
effective use of Lagrangian and Eulerian codes respectively. However conventional particle- 
in-cell and smooth particle hydrodynamics methods employ particles which are actually 
moving interpolation points. A new particle-based modeling methodology, termed 
Hamiltonian particle hydrodynamics, was developed by Fahrenthold and Koo (1997) to 
provide an alternative, fully Lagrangian, energy-based approach to shock physics 
simulations. This alternative formulation avoids the tensile and boundary instabilities 
associated with standard smooth particle hydrodynamics formulations and the diffusive grid- 
to-particle mapping schemes characteristic of particle-in-cell methods. 

In the work described herein, the method of Fahrenthold and Koo (1997) has been 
extended, by coupling the aforementioned hydrodynamic particle model to a hexahedral finite 
element based description of the continuum dynamics. The resulting continuum model 
retains all of the features (including general contact-impact effects) of Hamiltonian particle 
hydrodynamics, while in addition accounting for tensile strength, plasticity, and damage 
effects important in the simulation of hypervelocity impact on orbital debris shielding. A 
three dimensional, vectorized, and autotasked implementation of the extended particle 
method described here has been coded for application to orbital debris shielding design. 
Source code for the pre-processor (PREP), analysis code (EXOS), post-processor (POST), 
and rezoner (ZONE), have been delivered separately, along with a User's Guide describing 
installation and application of the software. 
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1. Introduction 

The Hamiltonian particle dynamics model developed by Fahrenthold and Koo 
(1997) is purely hydrodynamic. In order to account for strength effects important in the 
simulation of hypervelocity impacts on space structures, the aforementioned model has 
been extended, to incorporate elastic-plastic-damage effects. These effects are quantified 
using a three dimensional, hexahedral finite element characterization of deformation in the 
modeled continuum. 

The extension of the particle model of Fahrenthold and Koo (1997), described here, 
is based on a body centered cubic packing scheme for the material particles. Under such a 
packing scheme, each particle has eight nearest neighbors in the reference configuration. 
The center of mass coordinates of these eight nearest neighbors are designated as nodal 
coordinates for a hexahedral finite element centered initially on the particle. The volume of 
this element is calculated at each time step (using one point integration), and appears as a 
mechanical variable in the system level internal energy function. The relative velocity of the 
body centered particle, with respect to each of its nearest neighbors, is calculated at each 
time step to determine the local plastic strain rate. Finally a scalar continuum damage 
variable is introduced for each element, allowing for the loss of cohesion of the element 
when a user-specified level of plastic strain has been accumulated. 

The sections which follow describe in detail the hexahedral finite element based 
elastic-plastic-damage model implemented in the code EXOS, and illustrate application of 
the code by the simulation of example hypervelocity impact problems. A User's Guide for 
the analysis code EXOS, pre-processor PREP, post-processor POST, and rezoner ZONE 
is provided separately. In the interest of brevity the current report describes only the new 
finite element based augmentations developed to extend the model of Fahrenthold and 
Koo(1997). The reader is referred to the latter reference for a detailed description of the 
basic Hamiltonian particle hydrodynamics model. 



2. Kinematics 


The hexahedral element kinematics are similar to those employed in Lagrangian 
hydrocodes, for example DYNA3D (Hallquist, 1983). The components of the deformation 
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where c^ ^ denotes the 'k'th component of the center of mass position vector for nearest 
neighbor *j' of particle 'i'. The current (effective) volume for the element is then 

V eff(l) = 8 det[ f (l) ] (2) 


where 'det' denotes the determinant and the deformation gradient matrix is _f ^ . 
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3. Internal energy 

The total internal energy associated with particle 'i' is assumed to consist of particle 

and element based components, and take the form 

U = (1/2) I mW [u^CvW, sW) + (l -DW)uW(eW, sW ) ] (3a) 

i=l 

e (i) = v eff «/m(i) (3b) 

where m^, s^, D^, and are the mass, internal energy per unit mass, 

volume per unit mass, entropy per unit mass, continuum damage, and element volume 
associated with particle i'. Since the element volume depends on the particle coordinates 
c(i) and since the particle entropy (S^) and deformation gradient (F^) are related to the 
previously defined variables by 

SO) = m (i) s (i) ; v^/v® = det(F (l) ) = f(*) 3 (3c,d) 

where is the specific volume in the reference configuration, it follows that the system 

internal energy can be expressed as a function of the generalized coordinates F®, D^, 
and S® 


U = U(c®, F<3), sW, D(i)) 

(3e) 

Hence the generalized conservative forces for the system are 
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where iX*) is a damage energy release rate, is the particle pressure. 


p(j) du^X V^, s^ ) 
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(5a) 

and the particle and element temperatures are 
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The generalized conservative forces in equation (4a), associated with the 

particle center of mass coordinates, may be calculated from the element pressure (P e ^^), 
which is 


peff(i) = du (l) ( e (l) > ) 

de(») 

and the minor determinants of the element deformation gradient, which are 
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Equations (6) and (7) then determine the generalized forces for nearest neighbor n 1 (i) as 
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The generalized forces for nearest neighbor n2(i) are 
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The generalized forces for nearest neighbor n3(i) are 
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The generalized forces for nearest neighbor n4(i) are 
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The generalized forces for nearest neighbor n5(i) are 
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The generalized forces for nearest neighbor n6(i) are 
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The generalized forces for nearest neighbor n7(i) are 
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Finally the generalized forces for nearest neighbor n 1 (i) are 
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Note that when the damage for the element associated with particle 'i' is zero, the preceding 
generalized forces are zero, that is the element looses its cohesion. 

4. Plasticity model 

Just as the element's nodal displacements determine its current density and hence 
(in part) its internal energy, the relative velocities of the nodes with respect to the body 
centered particle determine the local rate of plastic strain. To be specific, the local plastic 
strain rate at particle 'i' is determined from the tangent relative velocities of particle 'i' with 
respect to its eight nearest neighbors in the reference configuration 

v tan(i,j) _ c(i) . c(j) - v^d) ; j = 1 , 2, ... 8 (11a) 

where 

v rad(i,j) _ { [ (c (i) . c (j)) / | c (i) _ c (j)| ] . (*0) - 6<J)) ) [ (c® - cO)) / Ic® - c0)| ] (lib) 

Given a yield stress (a^®) and current volume (V®) for the particle, the force on particle 
T associated with plastic flow is 



(11c) 


f P(‘> = X fPW) signlv^W] 
k j =1 k 

where 

f POJ) _ [ i / (h(i) + h(j))] [ i . minCD^), dO)) ] (1/2) ( v(0 <jy(0 + y(j) ay(0 ) (1 Id) 

The effective plastic strain (eP(0 ) is then calculated by integrating the rate equation 

g 

eP(0 = (2/3) Z { [ 1 / (h(0 + hO))] I v tan ('d) I ) 2 (lie) 

j=l 

for each particle. 

5. Damage evolution 

The damage evolution relation assumed here is the simplest possible, namely the 
damage is set to one when the plastic strain in an element reaches a user specified critical 
value. This value is normally termed an erosion or failure strain in the finite element 
literature. The introduction of more complex damage evolution relations and other failure 
criteria is relatively simple, and is under consideration for future versions of the analysis 
code. 

6. Entropy production 

The effects of plastic deformation and damage evolution must be accounted for in 
the entropy evolution relations, since plastic flow and damage evolution are dissipative 
processes. An irreversible entropy production rate for particle 'i* (s' rr ^) is calculated from 
the energy dissipation rate associated with plastic deformation and damage evolution, as 
follows: 

S irr(i) = (1/0(0 ) [ f P (l) . £(i) + rO) 6(0 ] (12) 

The coefficients of the generalized velocities in the entropy production relations will 
determine nonconservative generalized forces in the system level state equations. 
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7. Example simulations 

Appendices A, B, and C, show some example hypervelocity impact simulations 
performed using the code EXOS. 

The first example represents the oblique impact of a spherical projectile on a flat 
plate at 6.81 kilometers per second. Both the projectile and the shield were taken to be 
aluminum, with the material described by a Mie-Grtineisen equation of state. Parameters of 
the simulation are shown in Appendix A. Figures A-l, A-2, and A-3 show the simulation 
results. 

The second example represents the oblique impact of a cylindrical projectile on a 
flat plate at 7.0 kilometers per second. Again both the projectile and plate materials were 
taken to be aluminum, and a Mie-Grtineisen equation of state was used. Parameters of the 
simulation are shown in Appendix B. Figures B-l, B-2, B-3, and B-4 show the simulation 
results. 

The third example represents the oblique impact of a rod on a flat plate at 7.0 
kilometers per second. Again both the projectile and plate materials were taken to be 
aluminum, and a Mie-Grtineisen equation of state was used. Parameters of the simulation 
are shown in Appendix C. Figures C-l, C-2, C-3, and C-4 show the simulation results. 

8. Conclusion 

This report details the formulation and implementation of a finite element 
augmentation for the particle hydrodynamics model of Fahrenthold and Koo (1997). The 
coupling of particle based and finite element based models in a single code allows for the 
general characterization of both contact- impact effects and elastic-plastic-damage effects, 
important to the simulation of hypervelocity impacts on space structures. 



10 


References 


Fahrenthold, E.P., and Koo, J.C., "Hamiltonian Particle Hydrodynamics," 
Computer Methods in Applied Mechanics and Engineering, Vol. 146 (1997), pp. 43-52. 

Hallquist, J.O., 1983, Theoretical Manual for DYNA3D , Lawrence Livermore 
National Laboratory. 



APPENDIX A: Sphere impact simulation 


Simulation parameters: oblique impact of a sphere on a flat plate: 


Sphere diameter (aluminum) 

Impact velocity 

Impact obliquity 

Plate thickness (aluminum) 

Equation of state type 
Failure strain 
Yield stress 

Numerical shear viscosity coefficient 
Numerical bulk viscosity coefficient 
Numerical conduction coefficient 
Penalty stiffness coefficient 
Time step coefficient 
Number of particles 
Total simulation time 
Number of time steps 
CPU time (Cray J916) 


1.00 cm 

6.81 cm/psec 
15 degrees 
0.16 cm 
Mie-Griineisen 
0.50 

0.0029 Mbar 
0.01 
0.10 
1.00 

10.0 
10.0 
17,530 
3.00 |isec 
3,210 

2.82 hours 


List of figures (attached): 

Figure A-l. Oblique view: simulation a t = 0.0 psec. 
Figure A-2. Oblique view: simulation a t = 3.0 psec. 
Figure A-3. Normal view: simulation a t = 3.0 psec. 
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APPENDIX B: Cylinder impact simulation 


Simulation parameters: oblique impact of a cylinder on a flat plate: 


Cylinder diameter (aluminum) 

Cylinder length 

Impact velocity 

Impact obliquity 

Plate thickness (aluminum) 

Equation of state type 
Failure strain 
Yield stress 

Numerical shear viscosity coefficient 
Numerical bulk viscosity coefficient 
Numerical conduction coefficient 
Penalty stiffness coefficient 
Time step coefficient 
Number of particles 
Total simulation time 
Number of time steps 
CPU time (Cray J916) 


= 0.20 cm 

= 0.40 cm 

= 7.00 cm/psec 

= 3 1 degrees 

= 0.04 cm 

= Mie-Griineisen 
= 0.50 

= 0.0029 Mbar 

= 0.01 

= 0.10 

= 1.00 

= 10.0 

= 10.0 

= 11,184 

= 2.00 psec 

= 7,772 

= 4.66 hours 


List of figures (attached): 

Figure B-l. Oblique view: simulation a t = 0.00 psec. 
Figure B-2. Oblique view: simulation a t = 1.21 psec. 
Figure B-3. Oblique view: simulation a t = 2.00 psec. 
Figure B-4. Normal view: simulation a t = 2.00 psec. 
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APPENDIX C: Rod impact simulation 


Simulation parameters: oblique impact of a rod on a flat plate: 


Cylinder diameter (aluminum) 

Cylinder length 

Impact velocity 

Impact obliquity 

Plate thickness (aluminum) 

Equation of state type 
Failure strain 
Yield stress 

Numerical shear viscosity coefficient 
Numerical bulk viscosity coefficient 
Numerical conduction coefficient 
Penalty stiffness coefficient 
Time step coefficient 
Number of particles 
Total simulation time 
Number of time steps 
CPU time (Cray J916) 


0.08 cm 
0.24 

7.00 cm/psec 
31 degrees 
0.04 cm 
Mie-Griineisen 
0.50 

0.0029 Mbar 
0.01 
0.10 

1.00 
10.0 
10.0 
23,754 
1.00 psec 
7,959 
8.04 hours 


List of figures (attached): 

Figure C-l. Oblique view: simulation a t = 0.000 psec. 
Figure C-2. Oblique view: simulation a t = 0.206 psec. 
Figure C-3. Oblique view: simulation a t = 0.627 psec. 
Figure C-4. Oblique view: simulation a t = 1.000 psec. 
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